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Abstract 

Atmospheric conditions at the site of a cosmic ray observatory must be known for reconstructing 
observed extensive air showers. The Global Data Assimilation System (GDAS) is a global at- 
mospheric model predicated on meteorological measurements and numerical weather predictions. 
GDAS provides altitude-dependent profiles of the main state variables of the atmosphere like 
temperature, pressure, and humidity. The original data and their application to the air shower 
reconstruction of the Pierre Auger Observatory are described. By comparisons with radiosonde 
and weather station measurements obtained on-site in Malargiie and averaged monthly models, 
the utility of the GDAS data is shown. 

Keywords: Cosmic rays, extensive air showers, atmospheric monitoring, atmospheric models 
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1. Introduction 

The Pierre Auger Observatory [l|, |2J is located near the town of Malargiie in the province 
of Mendoza, Argentina. At the site, at the base of the Andes mountains, two well-established 
measurement techniques are combined to measure extensive air showers with energies above some 
10 17 eV. The hybrid detector consists of a Surface Detector (SD) array and five Fluorescence De- 
tector (FD) buildings. Each of the slightly more than 1600 SD stations is a water-filled Cherenkov 
detector, measuring the secondary particles of air showers that reach the ground. The detectors 
of the array are spaced by 1.5 km (750 m in a small infill area in the western part of the array) 
and provide the lateral particle distribution around a shower core. Four FD buildings comprise 
six telescopes each and one FD enhancement installation consists of three telescopes. In each FD 
telescope, the UV light emitted by excited nitrogen molecules along the shower track is collected 
by a large segmented mirror and reflected onto a camera composed of 440 PMTs. With this 
measurement, the geometry and the longitudinal profile of the shower can be obtained. 

For the reconstruction of extensive air showers, the optical properties of the atmosphere at 
the site of the observatory have to be known. This is particularly true for reconstructions based 
on data obtained with the fluorescence technique Q, but also impacts upon data collected with 
the surface detectors Q. The detection of clouds is an important task of the atmospheric mon- 
itoring systems. Clouds can obstruct or - through scattering of the intense Cherenkov light - 
amplify the apparent fluorescence light before it reaches the FD. To eliminate data recorded in 
cloudy conditions from physics analyses, lidar stations and infrared cloud cameras are installed 
at each FD station of the Pierre Auger Observatory. These instruments scan the fields of view 
of the fluorescence detectors several times per hour during data taking periods to measure the 
cloud coverage and the base height of clouds [5] . The vertical profile of the aerosol optical depth 
is measured once every hour using vertical laser shots from two facilities near the center of the 
array. Using the calibrated laser energy and the amount of light scattered out of the beam towards 
the FDs, the amount of aerosols can be estimated Weather conditions near ground, and the 
height-dependent atmospheric profiles of temperature, pressure and water vapor pressure are rele- 
vant for several Auger Observatory measurements. E.g., these parameters affect the production of 
fluorescence light by excited nitrogen molecules at the shower track, and the Rayleigh scattering of 
the light between the air shower and detector. Atmospheric conditions are measured by intermit- 
tent meteorological balloon radio soundings. Additionally, ground-based weather stations measure 
surface data continuously. The profiles from the weather balloons were averaged to obtain local 
models, called (new) Malargiie Monthly Models 0). Since March 2009, the atmospheric monitor- 
ing system has been upgraded with the implementation of a rapid monitoring system [6]. Part of 
the new program was the measurement of atmospheric profiles with radio soundings shortly after 
the detection of particularly high-energy air showers, a system called Balloon-the-Shower (BtS). 
This enables a high-quality reconstruction of the most interesting events. 

However, performing radio soundings and applying these data to air shower analyses is not 
straightforward. Very critical aspects are the time of the weather balloon ascent and the data 
validity period. Furthermore, performing radio soundings, in particular within BtS, imposes a 
large burden on the collaboration. Therefore, we investigate the possibility of using data from the 
Global Data Assimilation System (GDAS), a global atmospheric model, for the site of the Auger 
Observatory. The data are publicly available free of charge via READY (Real-time Environmental 
Applications and Display sYstem). Each data set contains all the main state variables with their 
dependence on altitude with a validity period of 180 minutes for each data set. 

Key aspects of the impact of the profiles of atmospheric state variables on the development 
and detection of extensive air showers are discussed briefly (Sec. [2]). We motivate the necessity 
of more reliable atmospheric profiles by a discussion about the data validity period of weather 
balloons (Sec. [3]), describe the content and processing of the GDAS data (Sec. |4j and compare 
them to local measurements (Sec. [5]). The new atmospheric data are implemented in the data 
processing and simulation framework of the Auger Observatory for an analysis of reconstructed 
air showers (Sec. [6]). 
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2. Impact of Atmospheric State Variables on the Development and Detection of Ex- 
tensive Air Showers 



Varying atmospheric conditions in terms of state variables like temperature, pressure and 
humidity, may alter the development and, in particular, the detection of extensive air showers. 
Here, different aspects relevant to the analysis of air showers at the Pierre Auger Observatory are 
discussed. 

The air fluorescence emission excited by the passage of an air shower depends on pressure, 
temperature, and humidity [7]. The collisional de-excitation of excited nitrogen molecules by 
other molecules of the atmosphere like nitrogen, oxygen, and water vapor counteracts the de- 
excitation of the molecules via radiation. These quenching processes are pressure and temperature 
dependent as described by kinetic gas theory, and dependent on the water vapor content in air. 
Furthermore, the collisional cross sections for nitrogen-nitrogen and nitrogen-oxygen collisions 
follows a power law in temperature, <j cx T a . Most recent experimental data indicate a negative 
exponent a. In reconstructions of air shower data from the Auger Observatory, the fluorescence 
yield with its dependence on atmospheric conditions is described using experimental results from 
the AIRFLY experiment [H, The absolute calibration of the main fluorescence emission at 
337.1 nm is taken from Nagano et al. [10j. The dependence of the fluorescence yield on atmospheric 
conditions translates to an atmospheric dependence of the reconstructed cosmic ray energy and 
the depth of shower maximum, the latter being an indicator for the mass of the primary cosmic 
ray particle. Even short-term variations of the atmosphere may introduce noticeable effects on 
these reconstructed parameters. 

Besides the fluorescence emission, the pressure, temperature and humidity profiles of the at- 
mosphere are important for other aspects of the reconstruction of data collected by the FD. These 
include the conversion between geometrical altitudes and atmospheric depth; the treatment of 
Cherenkov emission from air showers; and the transmission of the produced photons from the air 
shower to the FD: 

• The air shower development is governed by the interactions and decays of the secondary 
particles. These processes are largely determined by the atmospheric depth X, the total 
column density of atmospheric matter traversed by the air shower at a given point. X 
is calculated by integrating the density of air from the top of the atmosphere, along the 
trajectory of the shower through the gas. The observation of the longitudinal shower profiles 
by fluorescence telescopes is based on geometrical altitudes h. Thus, geometrical altitudes 
must be converted into atmospheric depth by taking into account the actual air density 
profile p(h) at the site of the Observatory, and the zenith angle 9 of the trajectory of the 
shower, 

cos J h 

• The secondary particles in extensive air showers travel faster than the speed of light in air. 
As a result, they induce the emission of Cherenkov light in a narrow, forward-beamed cone. 
Some of this light in the UV range may be - depending on the shower geometry relative to 
the FD telescope - detected together with the fluorescence light. To effectively subtract the 
Cherenkov photons from the total number of photons detected, the amount of Cherenkov 
light emitted by the air shower must be estimated. The Cherenkov yield depends on the 
refractive index n of the air, which itself depends on the wavelength of the emitted light 



as well as the temperature, pressure and humidity [llj, ll2( . Parameterized formulae for the 
refractive index of dry air, CO2 and water vapor are used to calculate a total refractive 
index, 

"tot - 1 = ("dry - 1) ' + (ncQ 2 - 1) • + («w ~ 1) ' — ■ (2) 

Pair Pair Pair 

The refractive index of each component is weighted with its density, which can be calculated 
using the number density and the molar mass of the constituent. Finally, the effect of the 



G 



decreasing number density with altitude is parameterized 13] as a function of pressure p 
and temperature $ in °C, 



l+p-(61.3-^).10' 10 

Hair L-[U tot 1J • f> • ggQgg 4 . Q _|_ Q. 003661 • 1?) ' W 

Between the production of fluorescence and Cherenkov light in the air shower and the de- 
tection at the FD telescope, the light is scattered by molecules in the atmosphere. The 



transmission of light depends on the Rayleigh cross section [14j , 



24-7T 3 /T7 2 — 1 > 

M\P,T,e) = jr - ] ^-(j^) •F aiz {\,p,e). (4) 



air 



where A is the wavelength in m and N the atmospheric molecular density, measured in 
molecules per m~ 3 . is the King factor that accounts for the anisotropy in the scattering 
introduced by non-spherical scatter centers, which depends slightly on pressure and humidity. 
The refractive index n air depends on several atmospheric state variables, see Eq. [3J 

The last three itemized effects on the reconstruction of extensive air showers can be taken 
sufficiently into account by using a proper description of the atmospheric state, e.g., the local 
atmospheric monthly models derived from multi-year meteorological balloon radio soundings (see 
Sec.|3|). They affect the reconstruction results of air shower data, mainly primary energy and posi- 
tion of shower maximum, only by marginally broadening the uncertainties without any significant 
systematic shifts. 

However, in the case of the earlier discussed fluorescence emission process and its atmospheric 
variability, systematic alterations of the reconstruction results may be seen together with increased 
uncertainties, even for short-term variations of the atmospheric parameters. Finally, after this 
discussion on atmospheric influences on FD analysis, it should be noted that uncertainties in the 
surface detector signals introduced by varying atmospheric conditions close to the ground are well 
understood and quantified 



3. Validity of Radio Soundings 

Since August 2002, meteorological radio soundings have been performed above the Pierre Auger 
Observatory to measure altitude-dependent profiles of atmospheric variables, mainly pressure, 
temperature, and relative humidity. Regular measurements were done until December 2008 in 
order to collect data for all months. After applying selection criteria, 261 profiles from the middle 
of 2002 until the end of 2008 could be used to build the new Malargiie Monthly Models 0. Starting 
in March 2009, the radio soundings became part of the rapid atmospheric monitoring system 
known as the Balloon-the-Shower (BtS) program [H,[l5j]. A fast online air shower reconstruction 
with subsequent quality selections is used to trigger the launch of a weather balloon by a local 
technician. 

A procedure was developed to find the period of time for which the data measured during 
the ascent of a weather balloon give a good description of the atmospheric conditions at the 
Pierre Auger Observatory. The 3-dimensional atmospheric conditions before and after a weather 
balloon ascent are unknown but data from local weather stations may help to identify stable 
periods or trends towards rapidly changing conditions. Every active weather station is used as an 
independent source of data, no matter how many stations contribute information during the period 
of the weather balloon ascent. For each station, the maximum variations of the temperature, the 
pressure, and the humidity are obtained for the duration of the corresponding weather balloon 
launch defined as the time between the start of the weather balloon and the burst of the balloon, 
see Fig. [TJ 

The difference AQ - with Q being temperature T, pressure p, or water vapor pressure e - 
between maximum (Q max ) and minimum (Q m i n ) values of every station during weather balloon 
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Figure 1: Schematic drawing of the procedure to find an extended time period of validity for data from 
weather balloons based on data measured by ground-based weather stations for typical atmospheric con- 
ditions. For details see text. 



flights can be seen in Fig. [2j From these histograms, periods with very stable conditions (AQ < 
AQi ow ), with typical conditions, and with unstable conditions are denned for each quantity Q (see 
caption of Fig. [2]). 




AT (K) Ap (hPa) Ae (hPa) 



Figure 2: Differences between maximum and minimum values for temperature T, pressure p, and water 
vapor pressure e as found in weather station data during a weather balloon launch. For temperature 
data (left), typical conditions are defined as variances between 4 K and 9 K, with 4 K being AQi ow for 
the quantity temperature. Periods with AT < 4 K indicate very stable conditions while periods with 
AT > 9 K are unstable. For pressure data (middle), typical conditions are defined between 1 and 3 hPa. 
For water vapor pressure data (right), the typical conditions are between 0.8 and 4.0 hPa. The very stable 
and unstable conditions of the last two quantities are defined accordingly. 



For typical conditions, the data of the weather stations are scanned before and after the 
time of the balloon ascent for each quantity Q, and the time at which any quantity leaves the 
range between Q min and Q ma x is determined. This time period gives the validity time period of 
the radio sounding for every active weather station. For launches performed during very stable 
conditions, the differences in weather station data are quite small. Thus, only small variations 
beyond the narrow interval would indicate the end of validity, imposing very strict cuts on this 
type of launch. For unstable conditions, the large AQ values could result in quite long extended 
periods of validity. Since both cases result in inappropriate validity periods, two special criteria 
for each quantity are found in addition to the typical case. For very stable conditions, Q max /min 
are redefined to <5 ma x/min = Q ± Qiow/2, where Q is the mean of the interval Q m i n to Q max and 
Qi ow is 4 K in the case of temperature data, 1 hPa for pressure data, and 0.8 hPa for water vapor 
pressure data. After definition of Q m ax/mim the same procedure as for the typical conditions is 
applied. In case of unstable conditions, the validity time period is set to the time period during 
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which the weather balloon ascended. 

The average duration of a weather balloon ascent was about 100 minutes. A validity time 
period of 200 minutes on average is given by the local weather station data as described above. 
Applying this procedure, about half of the cosmic ray events which triggered the BtS program are 
observed at times not covered by the period of validity of the corresponding balloon launch. 

Until its termination at the end of 2010, many details of local atmospheric conditions could 
be studied with the BtS program. The obtained atmospheric profiles can be applied to improve 
the reconstruction of the most interesting, high-energy air showers. However, the data are not 
suitable for application to the standard reconstruction because of their short period of validity. 
Only very few air shower events would be covered by atmospheric profiles from radio soundings. 

4. Global Data Assimilation System (GDAS) 

In the field of Numerical Weather Prediction, data assimilation is the adjustment of the de- 
velopment within a model to the real behavior of the atmosphere as found in meteorological 
observations (lij . The atmospheric models describe the atmospheric state at a given time and 
position. Three steps are needed to perform a full data assimilation: 

1. Collect data from meteorological measuring instruments placed all over the world. These 
instruments include weather stations on land, ships, and airplanes as well as radiosondes 
and weather satellites. 

2. Use a short-term forecast from a previous iteration of the numerical weather prediction to- 
gether with the measurements to describe the current situation. This additional information 
is needed because the available observations alone are not sufficient. The forecast or first 
guess adds more information to the system, namely all knowledge of atmospheric behavior 
expressed in mathematical model equations. The models use non-linear differential equations 
based on thermodynamics and fluid dynamics. 

3. Adjust the model output to the measured atmospheric state. The resulting 3-dimensional 
image of the atmosphere is called analysis. 

A schematic showing the principle of data assimilation is given in Fig. [3] At a given time to, 
the observations provide the value of a state variable. A model forecast for this variable from a 
previous iteration exists for the same time. The analysis step combines observation and forecast 
to describe the current state better than the forecast. This analysis is the initial point for the 
weather prediction model to create the forecast for a later time t\. 

4.1. GDAS Data 

The Global Data Assimilation System [l|| is an atmospheric model developed at NOAA'13 
National Centers for Environmental Prediction (NCEP). It provides an analysis four times a day 
(0, 6, 12, and 18 UTC) and a 3-, 6- and 9-hour forecast. The numerical weather prediction model 
used in the GDAS is the Global Forecast System (GFS). 

3-hourly data are available at 23 constant pressure levels - from 1000 hPa (roughly sea level) 
to 20 hPa (rj 26 km) - on a global l°-spaced latitude-longitude grid (180° by 360°). Each data 
set is complemented by data for the surface level. The data are stored in weekly files and made 
available online [l8(. In Table Q] the level indices corresponding to each data level are listed. For 
reference, the altitude from the US Standard Atmosphere 1976 (USStdA) [lj| is also given in the 
table. The actual height of the pressure level is stored in the data file. GDAS data are available 
starting January 2005. There are two periods without data in the sets. The first two weeks of May 
2005 and weeks 3 and 4 of November 2005 are missing. Other than that, the record is complete 
up to the present time (end of November 2011). 



National Oceanic and Atmospheric Administration. 
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Because of the lateral homogeneity of the atmospheric variables across the Auger array 
only one location point is needed to describe the atmospheric conditions. In Fig. |U the available 
GDAS grid points are marked as red crosses on a map together with a map of the surface and 
fluorescence detectors of the Auger Observatory. The grid point at 35° south and 69° west was 
chosen, at the north-eastern edge of the surface detector array. The two points to the west of 
the array are in the foothills of the Andes mountains and therefore not suitable. The point to 
the south-east of the array is quite far away and with a surface height of 1685 m a.s.L, it is also 
too high. Nevertheless, the profiles at this point are very similar to those at the chosen point, on 
average differing by less than 1 °C in temperature and less than 0.3 hPa in water vapor pressure 
at all altitudes, confirming the homogeneity 

The height at which the surface data are given changes over the years for the selected grid 
point. Starting in January 2005, the surface altitude is 1831.29 m above sea level. On May 31, 
2005, the surface height changes to 1403.38 m, and on August 22, 2006 it goes down further to 
1328.68 m and stays within a few centimeters of this value until July 27, 2010, when it changes 
to 1404.65 m. In Fig. [5l the surface height provided by the GDAS data sets is shown between 
January 2005 and December 2010. For reference, the altitudes of the lowest SD tank (1331.05 m) 
and the highest and lowest FD buildings (1712.3 m and 1416.2 m) are also shown. The reasons for 
these changes are regular improvements of the models and calculations used to produce the GDAS 
profiles, or resolution changes in the meteorological model [20(. These changes can occur again in 



the future, so the surface height of the data has to be monitored for undesired changes [21| . 

For air shower analysis, only data above ground level in Malargiie are interesting. Therefore, 
we only use data from the surface and from pressure levels 6 and above. The data from beginning 
of January to the end of May 2005 have a surface height of around 1800 m. This is even above the 
height of the highest FD building at Coihueco. We decided not to attempt an extrapolation down 
to the actual ground level of around 1300 m and discard these data. Therefore, the first data set 
we use is from June 1, 2005 at 0:00 UTC. 

4-2. Preprocessing of Data 

For air shower analyses, several types of information are stored in databases such as the one 
describing the state variables of the atmosphere. It contains values for temperature, pressure, 
relative humidity, air density, and atmospheric depth at several altitude levels. The first three 
quantities and the altitude are directly available in the GDAS data. Air density and atmospheric 
depth must be calculated. The surface data contain ground height, pressure at the ground, and 
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Figure 4: Map of the region around the array of the Pierre Auger Observatory with geographical latitudes 
and longitudes marked at the boundary of the figure. The positions of the surface and fluorescence 
detectors are superimposed. The available GDAS grid points for that area are marked as red crosses. The 
coordinates of Malargiie are —35.48° (south) and —69.58° (west), the selected GDAS grid point is —35° 
(south) and —69° (west). 



Table 1: Vertical levels of pressure surfaces. The actual height for each pressure level is given in the data 
file. For reference, the USStdA height is also listed in this table. Pressure level contains the surface 
values. 
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Figure 5: The surface height of the GDAS data between January 2005 and December 2010. For reference, 
the red line marks the height of the lowest tank in the array, and the blue lines represent the height of 
the FD buildings at Los Leones and Coihueco. 



relative humidity 2 meters above ground. Two temperature values are given, one at the surface 
and one 2 meters above ground. We decided to use the latter since we use it together with the 
relative humidity, which is also given 2 meters above ground, to calculate water vapor pressure. 

In the GDAS data, the altitude is given in geopotential meters with respect to a geoid (mean 
sea level) . In the air shower analysis framework of the Auger Observatory, geometric heights with 
respect to the WGS-84 ellipsoid are used. To move from geoid to ellipsoid, a constant value of 
26 m that arises from the geographic location of the Auger Observatory must be added to the 
height values of the model. The second step is to convert from geopotential height h to geometric 
altitude z (both measured in m), 

h 2 

z{h, $) = (1 + 0.002644 • cos(2$)) • h + (1 + 0.0089 • cos(2$)) • (5) 

OZ4GUUU 

where $ is the geometric latitude fjjj ]. 

To calculate the air density, the relative humidity must be converted into water vapor pressure, 
the partial pressure of water in air in Pa. This conversion depends on air temperature. The 
following approximation of the empirical Magnus formula is used in these calculations: 

e = -H_ . 6.1070 • exp ( R15 'M , > 0°C, (6) 
100% \234.9 + tf) ' y ' 

u ( 21.88-1? 



-ioo% 6 - 1064 - ex n^T^J' * < °° C ' (7) 

where u is the relative humidity in % and i? is the temperature in °C. Now, the air density in kg 
m -3 can be calculated with 

P ■ ^air 

= ~~R~T ri (8) 
where p is the pressure in Pa, T is the temperature in Kelvin, M a j r is the molar mass of air in 
kg mol^ 1 and R is the universal gas constant 8.31451 J K" 1 mol -1 . Moist air can be separated 
into three components to calculate its molar mass: dry air, water vapor and carbon dioxide. The 
molar mass of moist air is the sum of the molar masses of the components, weighted with the 
volume percentage of that component, 

M air = M dry • ip dry + M w ■ ip w + M C o 2 ■ yco 2 ■ (9) 
The molar masses for dry air, water vapor and CO2 are 0.02897, 0.04401 and 0.01802 kg mol -1 , 
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respectively. The volume percentage of CO2 is taken as 385 ppmv, the percentage of water (p w is 
the partial pressure e of water vapor divided by the pressure p, and dry air makes up the rest, 
fdry = 1 ^ <y5w ^ ¥>co 2 • The atmospheric depth can be calculated by integrating the air density in 
kg m~ 3 downward along a vertical track in the atmosphere, starting at infinity 

X(h )= / pdh. (10) 

J ho 

Here, the atmospheric depth is calculated in kg m~ 2 . In the astroparticle community, the unit g 
cm~ 2 is more common and will be used in this paper. 

Since the GDAS data only go up to around 25 to 30 km, we have to approximate the atmo- 
spheric depth at the top of the data profile using X — p/g and integrate numerically from that 
height down to ground level, g is the gravitational acceleration with dependence on altitude and 
geographical latitude [23j|, 

g(h,$)=g ($) ~ (3.085462 • 10~ 4 + 2.27 • 10~ 7 • cos(2$)) • h 
+ (7.254 • 10~ n + 1.0 • 10" 13 • cos(2$)) • h 2 
- (1.517-10" 17 + 6-10" 20 -cos(2$))-/i 3 , (11) 

where g is in cm s~ 2 and go ($) is the acceleration at sea level with dependence on latitude, 

ffo($) = 980.6160 • (1 - 0.0026373 ■ cos(2$) + 0.0000059 ■ cos 2 (2$)). (12) 

The integration is done by interpolating the density every 200 meters and using the trapezoidal 
rule to approximate the integral. 

For the simulation and reconstruction of air showers, the description of the atmospheric pa- 
rameters should ideally range from ground level to the top of the atmosphere. GDAS provides 
data between about 1400 and 30000 m. All profiles are extended up to 100 000 m - the approxi- 
mate boundary to outer space - using the US Standard Atmosphere 1976 [3] which describes the 
conditions above 30 000 m reasonably well. Below 1400 m, pressure, atmospheric depth, density, 
and water vapor pressure are exponentially extrapolated down to 1000 m based on the lowest two 
data points. The temperature profile is extrapolated linearly. Both extensions are outside the 
field of view of all FD stations. 

5. GDAS vs. Local Measurements 

To validate the quality of the GDAS data and to verify its applicability for air shower recon- 
structions at the Auger Observatory, GDAS data are compared with local measurements - atmo- 
spheric soundings with weather balloons and ground-based weather stations. The new Malargiie 
Monthly Models (nMMM) are also shown in some comparisons as a reference since they were the 
standard profiles used in reconstructions until recently. 

5.1. GDAS vs. Soundings with Weather Balloons 

Local atmospheric soundings have been performed above the array of the observatory since 
2002, but not on a regular basis. In the beginning, several week-long campaigns of launches were 
performed. Then, a pre-determined schedule was used to coordinate launches during and between 
dark periods with FD data taking and finally, with the BtS program, soundings are triggered by 
particularly high-energy events during FD data taking only. Most of the balloons were launched 
from the Balloon Launching Station (BLS, see Fig. |4|). To provide a set of atmospheric data for 
every measured air shower, the soundings were averaged to form monthly mean profiles. The 
current version of these Malargiie Monthly Models (nMMM) were compiled in early 2009. The 
models on average describe the atmosphere reasonably well, but show considerable fluctuations 
when compared to the actual sounding data 0] ■ The uncertainties of the profile for each variable 
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are given by the standard error of the variation within each month together with the absolute un- 
certainties of the sensors measuring the corresponding quantity. For atmospheric depth profiles, a 
picccwise fitting procedure is performed to ensure a reliable application of these parameterizations 
to air shower simulation programs. An additional uncertainty is included which covers the quality 
of the fitting procedure. 

In Fig. |6l the GDAS data are compared with the measured radiosonde data between 2005 
and 2008. The comparison between the radiosonde data and the monthly mean profiles is also 
shown. The monthly mean profiles are averaged from the sounding data until end of 2008, thus, 
they fit local sounding data of this period very well, see red squares in Fig. |6l The error bars 
denote the RMS of the differences at each height. The wave-like shape in the difference graph 
for the atmospheric depth X is driven by the piecewise parameterization of X in the nMMM. 
Above around 5 km, the GDAS data deviate only marginally from the measured data (black 
dots in Fig. |6]). Closer to the ground, slightly larger differences become apparent. The measured 
temperature is consistently higher than the model temperature. However, this might be caused 
by two problems acting in the same direction. For our local radio soundings, the temperature 
sensor might be quite often not properly acclimatized to outside conditions but is affected by the 
inside temperature of the Balloon Launching Station. It takes some minutes to overcome this 
effect during which time the weather balloon is already launched. When compiling the nMMM, 
only temperature data at 1600 m a.s.l. and above were taken from the radio soundings. To 
extrapolate to lower altitudes, a fit to these sounding data combined with data from local weather 
stations at 1401 m, 1420 m, 1423 m, 1483 m, and 1719 m was performed. However, data from 
the weather stations might also be influenced by the direct surface conditions beneath, since they 
are not standardized meteorological stations. It is also possible that the GDAS model does not 
adequately describe the heated surface of the elevated plain of the Pampa Amarilla and tends 
to assume free-atmosphere conditions. The pressure data of both models are in good agreement. 
The water vapor pressure fits almost perfectly, although the model values close to ground are both 
too low compared to real data which can be traced back to the difficult handling of humidity in 
general. In particular, the pressure differences propagate into the atmospheric depth where we see 
deviations from the measured data on the same scale as for the monthly mean profiles. 

In the comparison displayed in Fig. only radiosonde data from 2009 and 2010 are used in 
order to illustrate the strength of the GDAS model data. The nMMM are completely independent 
of this set of radio soundings. The most obvious changes compared with Fig. [5] are the worse 
descriptions of actual radio soundings with monthly mean profiles. Even though the uncertainties 
become larger, the nMMM do not describe the conditions measured during the years 2009 and 
2010 very well. In contrast, the GDAS data can represent the local conditions much better and 
the intrinsic uncertainty is consistently small. 

In the following, the non-perfect agreement between GDAS data and local sounding data below 
approximate 5 km a.s.l. is investigated in more detail. This effect close to ground could be caused 
by the nearby Andes and their influence on the climate above the array. GDAS was developed 
for global atmospheric predictions and therefore could be inadequate for very local atmospheric 
conditions. Under normal circumstances, the wind carries weather balloons launched at the Pampa 
Amarilla north-east or east with a horizontal displacement of about 100 to 150 km for ascents 
up to about 20 km a.s.l. In some cases however, the balloon just ascends with only a small 
horizontal drift in any direction, whereas the opposite extreme cases are horizontal displacements 
of more than 200 km. For this study, 25 short and localized ascents and 18 launches with very 
long balloon paths were selected manually. The differences between GDAS data and the measured 
data for both groups are shown in Fig. [8] The pressure data from GDAS do not describe these 
extreme conditions at the Auger Observatory as well as those seen under normal circumstances. 
For extremely short soundings, the local measurements reveal a high pressure area close to ground. 
This typically quite local effect is not reproduced by GDAS but for higher altitudes, the model 
data fit very well again. High pressure zones often indicate a stable atmospheric layering, so 
conditions change only on long-term scales. This could cause the good description of temperature 
profiles by GDAS. For the other extreme case, the launches with very long balloon paths, the local 
pressure data indicate a low pressure area, accompanied by significant winds. The conditions are 
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Figure 6: Difference of measured radiosonde data and GDAS models (black dots) and monthly mean profiles (red 
squares) versus height for all sounding data between 2005 and 2008. 
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Figure 7: Difference of measured radiosonde data and GDAS models (black dots) and monthly mean profiles (red 
squares) versus height for all radiosondes performed in 2009 and 2010. 
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dominated by turbulences, indicated by short-term and small-scale temperature variations. Thus, 
the GDAS model data do not fit the local measurements well. Overall, it is rather difficult to say 
if the topology of the Andes is the source of uncertainty or if it is also an effect of quite extreme 
weather conditions which might be induced by the structure of the Andes in the vicinity of the 
Pampa Amarilla. Moreover, both groups of investigated ascents are extremes and do not describe 
the usual conditions. Less than 20% of the ascents launched at the Pierre Auger Observatory fall 
into either category. 
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Figure 8: Difference of measured radiosonde data and GDAS models for short (black dots) and long ascents (red 
squares) versus height. 



Finally, possible inconsistencies between local measurements and GDAS data close to the sur- 
face are investigated by using weather station data. While radiosondes may suffer from measure- 
ment uncertainties near the ground, the weather stations are specifically designed to continuously 
measure the ground values. 

5.2. GDAS vs. Ground Weather Stations 

Five weather stations continuously monitor atmospheric values close to the ground, at about 2 
to 4 m above surface level. Four are located at the FD stations, and one was set up near the center 
of the array at the Central Laser Facility (CLF) . To make sure that the GDAS data describe the 
conditions at the ground reasonably well, the values provided by the GDAS data set are compared 
to all available weather station data. The GDAS data are interpolated at the height of the weather 
station. 

In the histograms on the left of Fig. |Hl the differences between measured weather station data 
and GDAS model data are shown for the weather stations close to the CLF and FD Loma Amarilla. 
All data measured in 2009 were used. Temperature, pressure, and water vapor pressure agree very 
well. Details of the histograms are listed in Tab. [2j 
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Figure 9: Difference of measured weather station data and GDAS models for all data from 2009. The mean and 
variance values for all histograms shown can be found in Table [2] 

On the left, the differences GDAS model minus measured data in temperature T, pressure p, and water vapor 
pressure e are shown for the weather station at the CLF (solid red line) and the station at FD Loma Amarilla 
(dashed black line). 

On the right, the differences between the data of the individual stations are shown. The solid black line represents 
the difference between data measured with the Los Leones (LL) weather station and the data from the CLF weather 
station. The solid green line corresponds to Los Morados (LM) minus CLF, the dashed red line is for Loma Amarilla 
(LA) minus CLF and the dotted blue line indicates Coihueco (CO) data minus CLF data. 
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Table 2: Mean and variance values for the histograms shown in Fig. [9] Values for the pressure were not 
corrected for the height difference of the stations, but they are consistent with those differences. Applying 
the USStdA, the pressure difference due to the different altitudes of the weather stations are —2.0 hPa 
(LL), -8.6 hPa (LA), -2.3 hPa (LM), and -32.8 hPa (CO). 
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On the right side of Fig. [9j the differences of the data of the individual weather stations 
are shown. The CLF weather station is close to the middle of the array and was chosen as a 
reference. Values for the pressure were not corrected for the height difference of the stations, 
but they are consistent with the height differences of the stations. The mean and width of the 
CLF-LA distribution is very similar to the distributions on the left for the GDAS data. These two 
histograms are expected to be similar because of the vicinity of the selected GDAS grid point and 
the FD building Loma Amarilla. Overall, the differences between the GDAS data and the weather 
station data are of the same order as the difference in data of two different weather stations. Only 
the difficult predictability of water vapor in the atmosphere close to ground can be seen again. 

In Fig. [10] the differences in temperature for the weather station at Los Morados are analyzed 
in more detail. At the top left, all data measured every 5 minutes are compared with GDAS data. 
Additionally, only station data every three hours at which new GDAS data are available are shown, 
scaled by a factor of 36. Both distributions show identical means of 0.1 K, as well as similar RMS 
values of 2.7 K and 2.6 K for all data and for the 3-hourly data, respectively. The atmospheric 
parameters at the Observatory are very stable during 3-hour periods. At the top right, the data 
are split into day and night. Night is defined as the UTC hours between and 10. Therefore, with 
14 compared to 10 hours of data, the daytime distribution contains more data. For Los Morados, 
both distributions are close to each other, with a mean of 0.0 K and RMS of 2.9 K during the day 
and a mean of 0.2 K and RMS of 2.4 K at night. A small difference at night is noticeable, with the 
GDAS data giving higher temperatures than the weather station. The distributions for different 
years (bottom left) and different seasons (bottom right), represented by one month of summer and 
one from winter, show no distinct features or differences. Similar studies for pressure and water 
vapor pressure mostly yield similar results. Nevertheless, two deviations from this general pattern 
are discussed in the following. 

In the left panel of Fig. [1TJ the temperature distributions for the station at Loma Amarilla 
are shown, split into daytime (solid red) and nighttime (dashed black) measurements. A clear 
separation of the two distributions is found. The mean of the distributions is —1.6 K for daytime 
and 1.6 K for nighttime, resulting in a difference of 3.2 K. This difference might be due to the local 
environment in which the weather station is placed. While some stations are far away from the FD 
buildings, others had to be mounted closer to or on top of other facilities. Thus, a standardized 
meteorological measurement can not be guaranteed. 

On the right of Fig. 111! the water vapor pressure for different seasons measured at the weather 
station at the CLF is shown. Clear differences are apparent. In austral winter (July), the water 
vapor pressure is very low, and the measurements agree better with GDAS data, the mean of the 
distribution is 0.8 hPa with an RMS of 1.2 hPa. In summer (January), the deviations are largest, 
where the mean difference drops to —1.4 hPa, and the RMS doubles to 2.4 hPa. This indicates 
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Figure 10: The difference of GDAS minus weather station data in temperature is shown for the weather station 
at Los Morados. At the top left, the difference of all data measured every 5 minutes (dashed black) is compared 
to the difference only of data measured every three hours when new GDAS data are available (solid red, scaled by 
36). At the top right, all data are split into daytime (dashed black) and nighttime (solid red) measurements. At 
the bottom left panel, all data are separated by years (2007 in solid blue, 2008 in solid red, 2010 in dashed black). 
At the bottom right, data are split into summer and winter, represented by one month each (January in dashed 
black, July in solid red). 




Figure 11: The difference of GDAS minus weather station data in temperature is shown in the left panel for the 
station at Loma Amarilla, separated by time of day (nighttime in dashed black, daytime in solid red). In the right 
histogram, the water vapor pressure for the station at the CLF is shown separated by seasons (Summer in dashed 
black, winter in solid red). 
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that the GDAS description of humidity is not perfect, underestimating the humidity in summer. 
However, the water vapor pressure calculation strongly depends on the temperature, so differences 
in temperature due to local effects of the surroundings of the station also affect this comparison. 

Apart from the differences that were seen between radiosonde data and GDAS data near the 
surface in the previous section, the comparison of GDAS data with weather station data shows 
a very good agreement. We conclude that the GDAS data describe the conditions at the Pierre 
Auger Observatory very well. Because of their highly reliable availability and high frequency of 
data sets, GDAS data are a suitable replacement for local radio soundings and also for the local 
monthly models. 

5.3. GDAS vs. Radio Soundings at Other Locations 

In addition to the measurements in Argentina, similar radiosonde launches were performed in 
south-east Colorado, United States. 27 weather balloons were launched between September 2009 
and December 2010 using identical radiosondes and equipment. The comparisons with the GDAS 
data for this location show differences of the same order as for the location of the Pierre Auger 
Observatory This is remarkable because most of the global atmospheric models, in particular 
those developed in Europe and North America, typically describe the conditions at the northern 
hemisphere much better. This is due to the fact that atmospheric measurements in South America 
and in general at the southern hemisphere are sparse and accurate modeling of the atmosphere is 
predicated on real data. 

Comparisons of GDAS data and radio soundings at further locations would go beyond the 
scope of this paper. Moreover, the radio soundings from south-east Colorado are independent 
data while other available data are from radio sounding databases which are part of the global 
meteorological network used for the creation of GDAS. 

6. Air Shower Reconstruction 

To study the effects caused by using the GDAS data, all air shower data from the Auger Obser- 
vatory collected between June 1, 2005 and end of 2010 are used in a reconstruction analysis using 
the Offline software framework of the Pierre Auger Observatory [24j. The change of atmosphere 
description will mainly affect the reconstruction of the fluorescence data, c.f. Sec. [21 Therefore, we 
concentrate on this part in the following. It is known that varying atmospheric conditions alter 
the fluorescence light production and transmission [3j]. The transmittance of the actual atmo- 
sphere is regularly measured during FD shifts and made available for air shower reconstructions 
via databases. The light production has to be calculated analytically during the reconstruction 
procedure. Its strong atmosphere-dependence as described in Sec. [5] is applied in the air shower 
reconstruction analysis. 

6.1. Data Reconstruction 

The following analysis is based on three sets of reconstructions. The first set, FY, is the 
until recently standard reconstruction of the Pierre Auger Observatory. The fluorescence yield 
is calculated with its atmosphere-dependence as described in 25], along with the monthly mean 
profiles (nMMM) obtained for the site of the Auger Observatory. For the second set, FY mo d, 
all currently known atmospheric effects in the fluorescence calculation are taken into account. 
Together with the standard atmosphere-dependence, the temperature-dependent collisional cross 
sections and humidity quenching are included (Q and references therein). Parameterizations for 
these two effects are taken from AIRFLY Q and later conference contributions by the AIRFLY 
collaboration. Again, the nMMM are used. The third set, FY™^ S , also explores the efficiency 
of the full atmosphere-dependent fluorescence description, but here the atmospheric nMMM are 
exchanged with the new 3-hourly GDAS data. 

Comparing the reconstruction sets with each other, the variation of the reconstructed primary 
energy E of air showers and the position of shower maximum A max can be determined, see Fig. 1121 
In the two upper figures, the binned difference of E and A max is displayed, and the dependences 
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on energy and month of these differences are shown in the figures in the middle and bottom, 
respectively. 

Using GDAS data for the reconstruction instead of nMMM affects the reconstructed primary 
energy only slightly. The mean of the difference FY™^ S minus FY mod is 0.4% with an RMS of 1.4% 
(Fig. [T21 top left, solid black line). For the reconstructed X max , only a small shift of —1.1 g cm~ 2 
is found with an RMS of 6.0 g cm~ 2 (Fig. [121 top right, solid black line). However, comparing 
the full atmosphere-dependent reconstruction FY™^ S with the unmodified reconstruction FY, a 
clear shift in E can be seen. An increase in E by 5.2% (RMS 1.5%) and a decrease of X max by 
— 1.9 g cm~ 2 (RMS 6.3 g cm~ 2 ) is found. These modified fluorescence settings are now used in 
the reconstruction of the Pierre Auger Observatory, in conjunction with other improvements to 
the procedure, see (26j |. 

The difference in reconstructed E vs. mean E (Fig. [12] middle left) reveals a small energy 
dependence, increasing towards higher energies. The description of atmospheric conditions close 
to the ground is very difficult in monthly mean profiles since the fluctuations in temperature 
and humidity are larger in the lower levels of the atmosphere (below 4 km) than in the upper 
layers. Consequently, a more precise description of actual atmospheric conditions with GDAS 
will alter the energy reconstruction compared with nMMM-based reconstructions for air showers 
that penetrate deeply into the atmosphere, usually high-energy events. The full atmosphere- 
dependent fluorescence calculation alters the light yield for conditions with very low temperatures, 
corresponding to higher altitudes. The energy dependence of the X max differences is a combined 
effect of slightly changed humidity conditions close to ground and temperature conditions higher up 
in the atmosphere together with the full atmosphere-dependent fluorescence calculation (Fig. [T2l 
middle right). 

The difference in energy is quite uniform throughout the year, see Fig. Q21 bottom left. For 
switching on GDAS instead of nMMM (black dots), it is confirmed that GDAS describes the 
conditions at the Auger Observatory very well, as good as the nMMM. Switching on the full 
atmosphere-dependent fluorescence calculation does not show a monthly dependence because the 
overall integral of the longitudinal light profiles is hardly changed, see e.g. [271 ] ■ Only the modi- 
fication of the shape of the longitudinal light profile causes a small monthly dependence of A max 
(Fig. [121 bottom right). 

In the following, some systematics caused by the particular shower geometry are studied. In 
the first set of figures (Fig. Q2J top row), the difference in E and X max vs. zenith angle of the 
shower 9 is displayed. The energy variation is quite uniform around the mean value up to about 
60°. Only more inclined showers show a stronger shift of reconstructed energy for the modified 
fluorescence yield calculation. Concerning the position of shower maximum, a dependence on 
zenith angle can be seen above 30°. 

Shown in the bottom row of Fig. [T31 the E and A max dependence on geometrical height of 
shower maximum gives a more pronounced view of the atmospheric conditions in combination 
with the atmosphere-dependent fluorescence description. Showers reaching their maximum at an 
average altitude between 3 and 7 km show the mean E difference as expected from Fig.[l2j top left. 
However, showers with very shallow or very deep positions of shower maximum are reconstructed 
with a 7-8% higher primary energy compared with that using the standard fluorescence calculation. 



The reconstructed A max follows the expectations according to the study shown in 27j . 



6.2. Impact on Shower Reconstruction Uncertainties 

To study the effect that the GDAS data have on the uncertainties of air shower reconstructions, 
air showers induced by protons and iron nuclei are simulated using the CONEX shower genera- 
tor [28[ with the QGSJETII hadronic interaction model [29[ for shower energies between 10 17 5 eV 
and 10 20 eV. The fluorescence light is generated including temperature-dependent collisional cross 
sections and vapor quenching. The time stamps of the air shower events correspond to the times 
of 109 radio soundings between August 2002 and December 2008 so that actual atmospheric pro- 
files can be used in the simulation. All 109 launches were performed at night during cloud-free 
conditions. After the atmospheric transmission, the detector optics and electronics arc simulated. 
The resulting data are then reconstructed using the radiosonde data, as well as the GDAS data. 
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Figure 12: Difference of reconstructed E (left) and X max (right). Top row: binned differences. Middle row: 
differences as a function of energy. Bottom row: differences as a function of month. Black solid line or filled dots 
for FY^^ S minus FY mot j, and red dashed line or open squares for FYj^* s minus FY. 
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Figure 13: Energy difference (left) and X max difference (right) vs. zenith angle 8 of EAS (top) and vs. geometrical 
height of shower maximum (bottom). Black filled dots are for FYjjDAS minus FY mo( j, and red open squares for 



FY^DAS minus p Y . 



A basic set of quality cuts is applied. The shower maximum has to be in the observed part 
of the track, and the uncertainty in reconstructed energy and A max must be below 20% and 
40 g cm~ 2 , respectively. Also, the Gaisser-Hillas profile fit has to have a given quality, x 2 /Ndf < 2, 
and the fraction of Cherenkov light from the shower observed by the telescopes must be less 
than 50%. After applying all cuts, the values of energy and X max from the reconstructions with 
different atmospheric description are compared. The differences in these reconstructions yield the 
uncertainties that are introduced by the use of the GDAS data instead of the actual atmospheric 
profile. 

The same study has been performed to determine the uncertainties from the nMMM |30| . 
There, the systematic error is less than 1% in energy and less than 2 g cm' 2 in X max . In the 
energy range from 10 17 - 5 eV to 10 20 eV, random energy- dependent reconstruction uncertainties of 
±1% and ±5 g cm~ 2 for low energies up to ±2% and ±7 g cm" 2 for high energies were found. In 
the course of our new study, we compute the same uncertainties due to the nMMM again. The 
main difference between the current analysis and the previous one is the implementation of the 
temperature-dependent collisional cross sections and the humidity quenching in the calculation of 
the fluorescence yield. 
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Fi gure 14: Energy difference (left) and X max difference (right) vs. reconstructed FD energy for simulated showers. 
The influences of the monthly models (top) and of GDAS data (bottom) are shown. Error bars denote the RMS 
spread. 



In Fig. [Ml both results of the new study are shown, the uncertainties due to nMMM and those 
due to GDAS. A deviation from zero indicates a systematic error, and the error bars denote the 
RMS spread of all simulated events and are a measure of the reconstruction uncertainty due to 
this atmospheric parameterization. In the top panels, the influence of the nMMM compared to the 
real atmospheric parameters from radiosondes is shown. The results for the systematic shifts are 
in agreement with the previous study [30j |. with only the systematic shift for A max in the lowest 
energy bin being higher. The RMS spread in energy is ±0.9% for low energies and up to ±2.4% 
for high energies. In X max , the RMS is between ±4 g cm -2 for low energies and ±6.5 g cm -2 for 
high energies. 

The influence on the reconstruction due to GDAS data is shown in the bottom part of Fig. Q31 
The systematic shifts in energy are of the same order, below 1%, but of opposite sign. The shifts 
in A max are much smaller than for the nMMM, less than 0.5 g cm" 2 . The RMS spread is also 
considerably smaller, ±0.9% and ±2 g cm -2 for low energies, and ±1.3% and ±3.5 g cm" 2 for 
high energies. The energy uncertainty at low energies is comparable to the uncertainty introduced 
by the nMMM, but at high energies the uncertainty is reduced by almost 50%. For A max , the 
uncertainties in all energy bins are halved. 

This study of the reconstruction uncertainties further demonstrates the advantages of GDAS 
data over the monthly mean profiles. 
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7. Conclusion 



The reconstructions of air showers measured at the Pierre Auger Observatory have used a set 
of monthly mean profiles as the standard atmospheric description until recently. These profiles 
are averages from meteorological radio soundings performed at the site of the observatory over 
several years. The mean profiles describe the local conditions reasonably well, but cannot describe 
short-term variations in the atmosphere. Because of the large burden radio soundings impose on 
the collaboration, and their ambiguous duration of validity, data from the Global Data Assimila- 
tion System (GDAS) are a welcome substitute for atmospheric descriptions. GDAS data rely on 
established meteorological models and have an excellent time resolution of 3 hours. 

A direct comparison of GDAS data for the site of the Auger Observatory with local atmospheric 
measurements validates the adequate accuracy of GDAS data with respect to horizontal and ver- 
tical as well as temporal resolution. The suitable online publication of these data by NCEP allows 
for an easy and timely updating of atmospheric databases used at the Pierre Auger Observatory. 

With an air shower reconstruction analysis, the applicability of GDAS data to air shower 
analyses can be confirmed along with an improved accuracy with respect to atmospheric conditions. 
Also, the value of using an atmosphere-dependent fluorescence description has been demonstrated. 
Using simulated air showers, we show that the GDAS data significantly reduce the systematic 
errors and overall uncertainties in air shower reconstructions. 

Because of the results discussed in this study, the standard air shower analyses of the Pierre 
Auger Observatory are now applying atmospheric data from GDAS and the fluorescence descrip- 
tion FY™^ s as described in Sec. [6] 
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